Cardiac output is the product of heart rate (hr) and stroke volume (sv).
We will use a parametric bootstrap to make predictions of cardiac output as a function of log10(body mass) and habitat based on our random-effect models for hr and sv.
n_boot <- 1000
We will make sets of many (1000) predictions for each of the species that occur in both the hr and sv data sets.
# note NAMES HAVE TO BE SAME as passed to glmmTMB() when models were fitted (geesh)
mammal_hr <- readRDS('fitted-models/hr-data.RDS')
hr_mod <- readRDS('fitted-models/hr-re-model.RDS')
mammal_sv <- readRDS('fitted-models/sv-data.RDS')
sv_mod <- readRDS('fitted-models/sv-re-model.RDS')
cardiac_species <- intersect(pull(mammal_hr, animal), pull(mammal_sv, animal))
The two datasets have 24 species in common.
Make dataset for which to make predictions, containing these species.
cardiac_hr <- mammal_hr %>%
dplyr::filter(animal %in% cardiac_species) %>%
dplyr::select(animal, order, genus, species, log10.body.mass, log10hr, habitat) %>%
mutate(across(where(is.factor), factor))
cardiac_sv <- mammal_sv %>%
dplyr::filter(animal %in% cardiac_species) %>%
dplyr::select(animal, order, genus, species, log10.body.mass, log10sv, habitat)
For each row of the dataset, make a set of 1000 model predictions of cardiac output. These predictions will include uncertainty in the model parameters (intercept and slopes) as well as variation attributed to order/genus/species. (Residual variation is not included). The mass estimates are different in the different datasets; here we use the mass estimate from the dataset to which the model (hr or sv) was originally fitted.
First we define custom functions to make the predictions
predict_hr <- function(.){
predict(.,
newdata = cardiac_hr,
re.form = NULL,
type = 'response')
}
predict_sv <- function(.){
predict(.,
newdata = cardiac_sv,
re.form = NULL,
type = 'response')
}
Do the bootstrap(s)
if (!file.exists('fitted-models/hr_boot.RDS')){
hr_boot <- bootMer(hr_mod,
FUN = predict_hr,
nsim = n_boot,
re.form = NULL,
type = 'parametric'
)
saveRDS(hr_boot, 'fitted-models/hr_boot.RDS')
}else{
hr_boot <- readRDS('fitted-models/hr_boot.RDS')
}
if (!file.exists('fitted-models/sv_boot.RDS')){
sv_boot <- bootMer(sv_mod,
FUN = predict_sv,
nsim = n_boot,
re.form = NULL,
type = 'parametric')
saveRDS(sv_boot, 'fitted-models/sv_boot.RDS')
}else{
sv_boot <- readRDS('fitted-models/sv_boot.RDS')
}
Add boot results to datasets, combine to one dataset, and compute predicted cardiac output (and its median and 95% bootstrap CI):
cardiac_hr <- cardiac_hr %>%
mutate(boot_pred_hr = unlist(apply(t(hr_boot$t), 1, list), recursive = FALSE))
cardiac_sv <- cardiac_sv %>%
mutate(boot_pred_sv = unlist(apply(t(sv_boot$t), 1, list), recursive = FALSE))
cardiac_preds <- inner_join(cardiac_hr, cardiac_sv,
by = c('animal', 'order', 'genus', 'species', 'habitat')) %>%
mutate(boot_pred_cardiac = map2(boot_pred_hr, boot_pred_sv, function(.x, .y) mapply(prod, 10^.x, 10^.y))) %>%
mutate(boot_median_cardiac = map_dbl(boot_pred_cardiac, median, na.rm = TRUE),
boot_lo = map_dbl(boot_pred_cardiac, quantile, probs = 0.025, na.rm = TRUE),
boot_hi = map_dbl(boot_pred_cardiac, quantile, probs = 0.975, na.rm = TRUE),
mean_mass = ((10^log10.body.mass.x + 10^log10.body.mass.y) / 2) )
# units: ml/stroke * beats/min
plot results
cardiac_mass <- gf_point(boot_median_cardiac ~ mean_mass, data = cardiac_preds,
color = ~habitat,
text = ~animal) %>%
gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = cardiac_preds,
color = ~habitat) %>%
gf_labs(title = 'Expected Cardiac Output\nfor the species modeled',
y = 'Predicted Cardiac Output (mL/min)',
x = 'Body Mass (kg)') %>%
gf_theme(scale_color_manual(values = my_colors)) %>%
ggplotly(tooltip = 'text')
cardiac_mass
htmlwidgets::saveWidget(cardiac_mass, "figures/cardiac_mass.html")
cardiac_mass_log <- gf_point(boot_median_cardiac ~ mean_mass, data = cardiac_preds,
color = ~habitat,
text = ~animal) %>%
gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = cardiac_preds,
color = ~habitat) %>%
gf_refine(scale_x_continuous(trans='log10')) %>%
gf_refine(scale_y_continuous(trans='log10')) %>%
gf_labs(title = 'Expected Cardiac Output\nfor the species modeled',
y = 'Predicted Cardiac Output (mL/min)), Log10 scale',
x = 'Body Mass (kg), Log10 scale') %>%
gf_theme(scale_color_manual(values = my_colors)) %>%
ggplotly(tooltip = 'text')
cardiac_mass_log
htmlwidgets::saveWidget(cardiac_mass_log, "figures/cardiac_mass_log.html")
Notes: The figure could also, of course, be non-interactive.
We can also make predictions for a “random selection of hypothetical new species of various masses” instead if you want to get a graph of lines with symmetrical error bands.
Total ventilation is the product of breathing rate (fr) and tidal volume (vt).
We will use a parametric bootstrap to make predictions of total ventilation as a function of log10(body mass) and habitat based on our random-effect models for fr and vt.
We will make sets of many (1000) predictions for each of the species that occur in both the fr and vt data sets.
# note NAMES HAVE TO BE SAME as passed to glmmTMB() when models were fitted (geesh)
mammal_fr <- readRDS('fitted-models/fr-data.RDS')
fr_mod <- readRDS('fitted-models/fr-re-model.RDS')
mammal_vt <- readRDS('fitted-models/vt-data.RDS')
vt_mod <- readRDS('fitted-models/vt-re-model.RDS')
vent_species <- intersect(pull(mammal_fr, animal), pull(mammal_vt, animal))
The two datasets have 31 species in common.
Make datasets for which to make predictions, containing these species.
vent_fr <- mammal_fr %>%
dplyr::filter(animal %in% vent_species) %>%
dplyr::select(animal, order, genus, species, log10.body.mass, log10fr, habitat) %>%
mutate(across(where(is.factor), factor))
vent_vt <- mammal_vt %>%
dplyr::filter(animal %in% vent_species) %>%
dplyr::select(animal, order, genus, species, log10.body.mass, log10vt, habitat)
For each row of the dataset, make a set of 1000 model predictions of ventilation. These predictions will include uncertainty in the model parameters (intercept and slopes) as well as variation attributed to order/genus/species. (Residual variation is not included). The mass estimates are different in the different datasets; here we use the mass estimate from the dataset to which the model (fr or vt) was originally fitted.
First we define custom functions to make the predictions
predict_fr <- function(.){
predict(.,
newdata = vent_fr,
re.form = NULL,
type = 'response')
}
predict_vt <- function(.){
predict(.,
newdata = vent_vt,
re.form = NULL,
type = 'response')
}
Do the bootstrap(s)
if (!file.exists('fitted-models/fr_boot.RDS')){
fr_boot <- bootMer(fr_mod,
FUN = predict_fr,
nsim = n_boot,
re.form = NULL,
type = 'parametric'
)
saveRDS(fr_boot, 'fitted-models/fr_boot.RDS')
}else{
fr_boot <- readRDS('fitted-models/fr_boot.RDS')
}
if (!file.exists('fitted-models/vt_boot.RDS')){
vt_boot <- bootMer(vt_mod,
FUN = predict_vt,
nsim = n_boot,
re.form = NULL,
type = 'parametric')
saveRDS(vt_boot, 'fitted-models/vt_boot.RDS')
}else{
vt_boot <- readRDS('fitted-models/vt_boot.RDS')
}
Add boot results to datasets, combine to one dataset, and compute predicted ventilation (and its median and 95% bootstrap CI):
vent_fr <- vent_fr %>%
mutate(boot_pred_fr = unlist(apply(t(fr_boot$t), 1, list), recursive = FALSE))
vent_vt <- vent_vt %>%
mutate(boot_pred_vt = unlist(apply(t(vt_boot$t), 1, list), recursive = FALSE))
vent_preds <- inner_join(vent_fr, vent_vt,
by = c('animal', 'order', 'genus', 'species', 'habitat')) %>%
mutate(boot_pred_vent = map2(boot_pred_fr, boot_pred_vt, function(.x, .y) mapply(prod, 10^.x, 10^.y))) %>%
mutate(boot_median_vent = map_dbl(boot_pred_vent, median, na.rm = TRUE),
boot_lo = map_dbl(boot_pred_vent, quantile, probs = 0.025, na.rm = TRUE),
boot_hi = map_dbl(boot_pred_vent, quantile, probs = 0.975, na.rm = TRUE),
mean_mass = ((10^log10.body.mass.x + 10^log10.body.mass.y) / 2) )
# units: breaths/minute (?) * mL / breath ---> mL/min
plot results
vent_mass <- gf_point(boot_median_vent ~ mean_mass, data = vent_preds,
color = ~habitat,
text = ~animal) %>%
gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = vent_preds,
color = ~habitat) %>%
gf_labs(title = 'Expected Ventilation\nfor the species modeled',
y = 'Predicted Ventilation (mL/min)',
x = 'Body Mass (kg)') %>%
gf_theme(scale_color_manual(values = my_colors)) %>%
ggplotly(tooltip = 'text')
vent_mass
htmlwidgets::saveWidget(vent_mass, "figures/vent_mass.html")
vent_mass_log <- gf_point(boot_median_vent ~ mean_mass, data = vent_preds,
color = ~habitat,
text = ~animal) %>%
gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = vent_preds,
color = ~habitat) %>%
gf_refine(scale_x_continuous(trans='log10')) %>%
gf_refine(scale_y_continuous(trans='log10')) %>%
gf_labs(title = 'Expected Ventilation\nfor the species modeled',
y = 'Predicted Ventilation (mL/min), Log10 scale',
x = 'Body Mass (kg), Log10 scale') %>%
gf_theme(scale_color_manual(values = my_colors)) %>%
ggplotly(tooltip = 'text')
vent_mass_log
htmlwidgets::saveWidget(vent_mass_log, "figures/vent_mass_log.html")
Notes: The figure could also, of course, be non-interactive.
We can also make predictions for a “random selection of hypothetical new species of various masses” instead if you want to get a graph of lines with symmetrical error bands.